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ABSTRACT 

Computational fluid dynamics (CFD) has proven to be 
an invaluable tool for the design and analysis of high- 
speed propulsion devices. Massively parallel comput- 
ing, together with the maturation of robust CFD codes, 
has made it possible to perform simulations of com- 
plete engine flowpaths. Steady-state Reynolds-Averaged 
Navier-Stokes simulations are now routinely used in the 
scramjet engine development cycle to determine optimal 
fuel injector arrangements, investigate trends noted dur- 
ing testing, and extract various measures of engine ef- 
ficiency. Unfortunately, the turbulence and combustion 
models used in these codes have not changed significantly 
over the past decade. Hence, the CFD practitioner must 
often rely heavily on existing measurements (at similar 
flow conditions) to calibrate model coefficients on a case- 
by-case basis. This paper provides an overview of the 
modeled equations typically employed by commercial- 
quality CFD codes for high-speed combustion applica- 
tions. Careful attention is given to the approximations 
employed for each of the unclosed terms in the averaged 
equation set. The salient features (and shortcomings) of 
common models used to close these terms are covered in 
detail, and several academic efforts aimed at addressing 
these shortcomings are discussed. 

INTRODUCTION 

Computational Fluid Dynamics (CFD) models typi- 
cally employed for compressible reacting internal flows 
have far less predictive capabilities than their counter- 
parts used for low-speed external flow applications. CFD 
models and experimental techniques applied to low-speed 
external flows have reached a level of maturity such that 
commercial aviation companies are now asking for drag 
coefficient estimates that have an uncertainty level as low 
as one count 1 or 0.0001. This translates to an uncer- 
tainty level of ±0.5% of the total drag for a typical air- 
craft at cruise conditions. While this level of accuracy 
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has not yet been realized, a statistical analysis of drag 
predictions presented at a recent Drag Reduction Work- 
shop 2 organized by the American Institute of Aeronau- 
tics and Astronautics showed that drag was predicted by 
CFD to within an uncertainty of ±43 drag counts. This 
value is comparable to (although larger than) the uncer- 
tainty of ±8 drag counts extracted from results based on 
wind tunnel tests. This level of expectation from CFD 
data is far removed from that felt by CFD practitioners 
of high-speed reacting internal flows. Results presented 
at a recent Joint Army/Navy/NASA/Air Force Work- 
shop on Turbulence and Kinetics Models for Scramjet 
Simulation included several examples where Reynolds- 
Averaged Navier-Stokes (RANS) models failed to even 
qualitatively mimic the fundamental flow physics present 
in these devices. 

Higher order modeling approaches, such as Large 
Eddy Simulation (LES), offer significant advantages that 
overcome many of the shortcomings associated with the 
statistical representation of single-point RANS closures. 
The LES approach for turbulence closure attempts to 
resolve the large-scale components of turbulence while 
modeling the smaller scales. Most of the transport of 
mass, momentum, and energy (on the order of 90%) is 
done by the large eddies, while the primary role of the 
small eddies is to dissipate these fluctuations. Hence, 
it is the large eddies that tend to interact directly with 
the mean flow. The resolution of the large scales implies 
that values chosen for modeled turbulent transport coeffi- 
cients (e.g. turbulent Prandtl and Schmidt numbers) will 
have less of an impact on the overall flowfield prediction. 
The smaller turbulent scales tend to be isotropic in na- 
ture and less dependent on boundary conditions and flow 
type than the larger scales. Thus, the modeling developed 
for small scales should be more generally applicable than 
models developed for the entire range of turbulent scales. 
Unfortunately, the computational costs of LES often pro- 
hibit its use as an engineering design tool for practical 
applications. This is particularly true for attached wall- 
bounded regions at modest to high Reynolds numbers. In 
fact, Spalart 3 has estimated that the application of tradi- 


1 

American Institute of Aeronautics and Astronautics 


This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States. 



tional LES to an airliner wing would require on the order 
of 10 20 floating point operations. This value is roughly 
one million times that of the largest RANS calculations 
attempted today. 

The immense costs involved with resolving even a 
fraction of the turbulence spectra forces the contin- 
ued reliance on single-point phenomenological models 
for the foreseeable future. Therefore, enhancements 
to Reynolds- Averaged Navier-Stokes methodologies will 
continue to be in high demand. This paper summarizes 
the current state-of-the-art modeling procedures used by 
engineers to model high-speed reacting flows. The typi- 
cal set of equations used for high-speed propulsion appli- 
cations are described along with all of the approximations 
(many of which are often taken for granted) required to 
close the averaged equation set. Some recent academic 
works meant to expand the applicability of the modeled 
equation set are also highlighted. 

GOVERNING EQUATIONS 


The equations that describe chemically reacting single- 
phase flows at conditions representative of most high- 
speed combustion applications are the Navier-Stokes 
equations coupled with ns— 1 species mass continuity 
equations (ns is the number of species considered). These 
partial differential equations can be written as follows: 



J t + 37 (p Hu j + 4j ~ r ij u i) = 0 ( lc ) 


d d 

lh ( pKm ) + Jx { pYrnU j + P^,) = w m (Id) 

where p is the density, u i is the velocity, P is the pressure, 
E is the total energy, H is the total enthalpy, T i; is the 
stress tensor, qj is the heat flux vector, and Y m , Vj, and vv m 
are the mass fraction, diffusion velocity, and production 
rate, respectively, of species “m”. 

The time-averaged equations are obtained by decom- 
posing each flow variable into a mean and fluctuating 
part. The following combination of conventional 
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is a common choice that minimizes the number of un- 
known correlations that appear. Substituting the decom- 
posed variables (Eq. 4) into Eqs. la - Id and averaging 
the result yields the desired time-averaged equation set: 

a7 + ^K) =0 (5a) 
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All terms on the right-hand-side of Eqs. 5b - 5d require 
modeling assumptions. 

Two unclosed terms arise in the time-averaged mo- 
mentum equation (Eq. 5b). The first term is the time- 
averaged molecular stress tensor, and the second term 
is the Reynolds stress tensor. The Reynolds stress ten- 
sor is the predominant term, and nearly all of the mod- 
eling effort devoted towards the closure of the momen- 
tum equations has centered around this term. The classes 
of models available for this term will be described later. 
The remaining term (T-) has historically been modeled 
by ignoring the effects of turbulent fluctuations on the 
molecular viscosity, p, and assuming that the conven- 
tional (a t ) and mass-weighted («,) average velocity are 
approximately equal. For a Newtonian fluid, these as- 
sumptions allow the average stress tensor to be approxi- 


mated as follows: 
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Direct Numerical Simulation (DNS) studies 4 ' 5 have sup- 
ported assumptions of this type, at least for perfect gases 
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under mildly compressible conditions. 

The time-averaged energy equation (Eq. 5c) introduces 
three additional correlations that require modeling (t j; 
has already been considered). The first new term is a 
molecular diffusion term that is well approximated for 
incompressible flows by the following expression: 


shown later. The second term is the dot product of the 
mean velocity with the Reynolds stress tensor. This term 
is closed based on the model chosen for the Reynolds 
stress tensor. The third term represents turbulent trans- 
port of the turbulent kinetic energy. The gradient diffu- 
sion approximation is typically used to model this term, 



where k is the turbulent kinetic energy, 

k = ^u"u" (8) 

For compressible flows, one typically assumes that this 
relationship remains valid. The time-averaged heat flux 
vector usually contains contributions from heat conduc- 
tion and an energy flux due to inter-species diffusion, i.e. 

ZS’j' ns 
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where ju, is the eddy viscosity and a k is a closure coeffi- 
cient defined by the chosen model for turbulence. 

The remaining terms that require closure reside in the 
species continuity equations (Eq. 5d). The first term is the 
species production rate. A multitude of models exist for 
closing this term that range from simple eddy break-up 
/ eddy dissipation models 6 ’ 7 to more elaborate methods 
based on probability theory. 8-10 A description of mod- 
els that are typically used for high-speed reacting flows 
will be discussed later. The remaining unclosed terms are 
diffusion terms. The diffusion velocity of species “m” is 
usually evaluated from Fick’s law of diffusion, i.e. 


The contribution from heat conduction is modeled in a 
manner that is consistent with that done for the molecu- 
lar stress tensor. That is, turbulent fluctuation effects are 
omitted when evaluating the thermal conductivity. A, and 
the mass-weighted temperature fluctuation average, T", 
is assumed to be negligibly small, i.e. 




(10) 


The treatment of the time-averaged inter-species diffu- 
sion term varies depending on the model chosen for the 
species diffusion velocity. The final term in Eq. 5c to be 
modeled is The average (mass-weighted) total 

enthalpy can be written in terms of the static enthalpy, h, 
and kinetic energy terms, 


H = h+ l -{u i u i + 2k) (11) 

Subtracting this expression from the expanded instanta- 
neous total enthalpy yields the fluctuating component of 
the total enthalpy, i.e. 

H H II _ 

H =h +u i u i +k (12) 

The unclosed correlation, pH" u’j, can then be expanded 
to yield 
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when the Reynolds-averaged equation set is considered. 
In this expression, D is the mass diffusivity of species 
“m” relative to the mixture. The use of this expression, 
in lieu of the costly evaluation of the multicomponent 
diffusion equation, is often justified by the premise that 
the “effective” turbulent diffusion is expected to domi- 
nate the molecular diffusion processes throughout most 
of the fiowfield. Through Fick’s law, the terms involving 
the species diffusion velocities are expressed as follows: 
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If one employs the same approximations used to model 
the average molecular stress tensor (i.e. turbulent fluctua- 
tion effects neglected on the mixture diffusivity, and con- 
ventional averages assumed equivalent to mass-weighted 
averages), then these expression simplify to the follow- 
ing: 
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The first term is the Reynolds heat flux vector. This term 
is modeled with various levels of complexity as will be 


Note that the effect of temperature fluctuations on the 
species enthalpy had to be ignored to arrive at the above 
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expression for the averaged inter-species diffusion terms. 
This rather dubious approximation (along with the ne- 
glection of composition fluctuations) is also made when 
extracting the mean temperature from the conserved vari- 
ables (or vice-versa). The final term to be modeled is the 
Reynolds mass flux vector, pT"u”. This term is usually 
modeled with the gradient diffusion approximation, al- 
though more complex models have also been used as will 
be shown later. 

The model chosen for the equation of state intro- 
duces additional closure uncertainties. Even the simplest 
choice, where the fluid is assumed to behave as a mixture 
of perfect gases, 

P = pRT (18) 

requires modeling assumptions since the gas “constant” 
( R ) varies with composition, 

ns y 
m= 1 

In this expression, R u is the universal gas constant and 
W m is the molecular weight of species “m”. All efforts, 
known to the author, circumvent this closure difficulty by 
simply neglecting the effects of composition fluctuations 
on the equation of state, i.e. 

P = ^T«pR(? m )f (20) 

MODELING PRACTICES 

The previous section highlighted the numerous closure 
approximations that are required when modeling com- 
pressible reacting flows. This section discusses common 
closure approximations to the terms that are generally 
perceived to require the most attention by model develop- 
ers. These terms are the Reynolds stress tensor (pu"u'j), 
Reynolds heat flux vector (phi' u”), Reynolds mass flux 

vector (pY^u'j), and the time-averaged chemical source 
term (vv m ). A description of models that are typically em- 
ployed in high-speed combustion applications and their 
known deficiencies are described in the sections that fol- 
low. 

Reynolds Stress Tensor 

The most common closures used for the Reynolds 
stress tensor are linear models based on the Boussinesq 
approximation, i.e. 


( ~ du t \ 

(da: 

duj \ 

r + *si) 


+ d7 i ) 


These models assume that the Reynolds stress compo- 
nents are related to the mean strain rate tensor through an 


isotropic eddy viscosity (p,). Models for the eddy vis- 
cosity vary in complexity from simple algebraic (zero- 
equation) models 1 1 which require specification of a tur- 
bulent velocity and length scale, to two-equation mod- 
els 12-15 which solve partial differential equations for 
both the turbulent velocity scale and an additional tur- 
bulence scale (e.g. a length scale, time scale, or dissi- 
pation rate). A three-equation ( k-e-v 2 ) model has also 
been proposed in the literature, 16 although this model 
has not been extensively applied to high-speed reacting 
flows. Algebraic models have the advantage of being nu- 
merically robust and easy to implement (at least for rel- 
atively simple geometries). However, these models of- 
ten require changes in their coefficients when applied to 
different types of flowfields, and ambiguities often arise 
when defining the turbulence scales for complex geome- 
tries. Two-equation models, on the other hand, tend to 
have a larger range of applicability, and they are easily 
extended to complex geometries where it may be diffi- 
cult to define relevant turbulent scales algebraically. One- 
equation models, that involve a transport equation for a 
quantity that can be directly related to the eddy viscos- 
ity, 17, 18 have gained popularity in recent years, particu- 
larly for external flow applications. This trend has not yet 
been seen, to a large degree, for internal reacting flows. 

The linear eddy viscosity models described previously 
have several deficiencies that are rectified by invoking 
higher order models. The first deficiency is a result of 
the direct proportionality assumed between the Reynolds 
stress and mean strain rate tensors (i.e. the Boussinesq 
approximation). This feature prevents the prediction of 
secondary flow motions that result from Reynolds stress 
anisotropies. Moreover, these models do not incorpo- 
rate the influences of pressure-strain correlations, which 
are responsible for the distribution of anisotropy among 
the normal stress components. The linear eddy viscosity 
models are also unable to rigorously account for stream- 
line curvature effects, since the Reynolds stresses depend 
solely on the frame-invariant strain rate tensor. These de- 
ficiencies are resolved through the use of second order 
models that involve transport equations for each of the 
Reynolds stress components, i.e. 





(V) 
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This expression shows that the Reynolds stress ten- 
sor is influenced by: (I) advection, (II) turbulent convec- 
tion, (III) source/sink due to the mean velocity gradient, 
(IV) source/sink due to the pressure gradient, and (V) dis- 
sipation due to viscosity. The unclosed terms are terms 
(II), (IV), and (V). The chain rule has historically been 
applied to term (IV), resulting in an additional diffusion 
term and a pressure-strain rate correlation. 


stress equation. Implicit algebraic models 25 are ob- 
tained by enforcing equilibrium assumptions on the tur- 
bulence. The specific assumptions are that the turbulence 
has reached an equilibrium state, i.e. 
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Similar expansions are applied to term (V) to yield a dif- 
fusion term and a dissipation term, i.e. 


„ dr jk , „dx ik 

dx k +U Jdx k 
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The pressure-strain rate correlation is responsible for the 
distribution of anisotropy among the normal stress com- 
ponents. This term is often of the same order of magni- 
tude as the source term due to the mean velocity gradient, 
term (III). Hence, substantial efforts have been devoted 
towards the modeling of this term. This term is usually 
partitioned into a “slow” relaxation towards isotropy term 
and a “rapid” response term resulting from imposed mean 
velocity gradients. 19-21 The dissipation term in Eq. 24 is 
typically partitioned into isotropic and deviatoric compo- 
nents, with the deviatoric component neglected in most 
works, i.e. 


du" du"j 2 ~ 


where 


£= pdu i du L 
P dx k dx k 


(25) 


(26) 


Models for the third order velocity correlation, term (II), 
can be found in Refs. 22 - 24. This term accounts for the 
turbulent convection of the Reynolds stress and is mod- 
eled to mimic a diffusion process. 

The computational cost associated with solving the 
Reynolds stress transport equations has discouraged its 
use for complex engineering calculations. The increased 
cost is due to the additional transport equations, and the 
stiffness posed by the highly non-linear relationships in- 
troduced to close these equations. This has led many to 
consider algebraic closures derived from the Reynolds 


and any anisotropies resulting from the turbulent trans- 
port and diffusion terms are proportional to anisotropies 
in the Reynolds stresses, 


dx. 


[P u "i u "j u "k + Pu "i 8 jk + Pu "j 5 ik ~ u "i x jk ~ u "j T ,k) = 
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Unfortunately, the iterative process required to solve the 
resulting system of equations has proven to be extremely 
“stiff”. Pope 26 was able to cast these expressions into a 
set of equations that result in explicit relationships for the 
Reynolds stresses for two-dimensional flows. This ap- 
proach was subsequently extended to three-dimensional 
flows. 27, 28 Algebraic Reynolds stress models, in contrast 
to linear eddy viscosity models, retain the information 
from the pressure-strain correlation models of the full 
Reynolds stress closure, and allow for Reynolds stress 
anisotropies. Of course, the applicability of the formu- 
lation hinges on the validity of the equilibrium assump- 
tions given by Eqs. 27 and 28. The computational ex- 
pense associated with explicit algebraic Reynolds stress 
models (EASM) is only slightly greater than that required 
for standard two-equation variants of linear eddy viscos- 
ity closures. 

The importance of accounting for Reynolds stress 
anisotropies can be illustrated by considering flow 
through a rectangular duct. These flowfields are known 
to contain stress-induced secondary motions near the cor- 
ners of the duct, which develop due to Reynolds stress 
anisotropies. Computed pitot pressure distributions ex- 
tracted from simulations of a Mach 3.9 flow in a square 
duct are compared with measurements 29 in Figs, and . In 
these figures, the measured data is shown on the left of the 
symmetry plane, while the computed results are shown on 
the right. Fig. compares computed results using the lin- 
ear Wilcox k -(0 model 13 to measurements, while Fig. 
compares computed results using an explicit algebraic 
Reynolds stress model 30 to measurements. When the 
Reynolds stress anisotropies are not accounted for (Fig. ), 
secondary flow structures do not develop in the comer re- 
gion of the duct. As a result, the boundary layer builds up 
more rapidly near the comers. The algebraic Reynolds 
stress model accounts for the stress anisotropies, allow- 
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ing the secondary flow structures to develop. These struc- 
tures transport high momentum fluid from the core flow 
into the corner regions which results in a “squared-off” 
boundary layer profile that more closely matches the ex- 
perimental data. The importance of accounting for the 
secondary flow structures is further illustrated in the wall 
pressure trace along the duct (see Fig. ). The boundary 
layer growth is over-predicted in the linear k -(0 results, 
yielding a larger total pressure loss in the latter half of 
the duct. Rectangular flowpaths are prevalent in many 
scramjet propulsion systems, hence accounting for the 
stress anisotropies may prove to be a critical ingredient 
when assessing inlet and isolator performance. 

As a final note, a few statements should be made for 
the class of models known simply as non-linear eddy vis- 
cosity models (NLEVM). 31-33 These models are func- 
tionally similar to EASM models, in the sense that the 
Reynolds stress tensor is represented by a polynomial ex- 
pansion of some given tensorial basis. 34 The primary dif- 
ference between EASM and NLEVM models is the man- 
ner in which the expansion coefficients are determined. 
The expansion coefficients derived for EASM models are 
based on the mathematical procedure followed to recast 
the implicit Reynolds stress expressions into explicit re- 
lationships. The expansion coefficients determined for 
NLEVM models, on the other hand, are based on empiri- 
cism and realizability constraints. 

Reynolds Heat/Mass Flux Vector 

The turbulent transport of a scalar property has histor- 
ically been modeled using the gradient diffusion hypoth- 
esis. This model choice assumes that the turbulent trans- 
port of the scalar is in the direction of decreasing value 
for that scalar. This leads to the following model expres- 
sions for the Reynolds heat flux 



P, dh 

Pr, d Xj 


and mass flux 


pY m “j 


Pt d? m 
Sc, dxj 


(29) 


(30) 


vectors. The diffusion rates are controlled by specifying 
the turbulent Prandtl (Pr,) and Schmidt (Sc,) numbers. 
The turbulent Prandtl number specifies the ratio of the 
rate of turbulent momentum transport to rate of turbu- 
lent energy transport, while the turbulent Schmidt num- 
ber defines the ratio of the turbulent momentum transport 
rate to turbulent mass transport rate. Constant values for 
these coefficients are usually assumed in applications for 
low- and high-speed reacting flows of engineering inter- 
est, even though values for these coefficients have been 
shown to vary spatially. 35-45 Table 1 summarizes the 


range of values that have been observed (both experimen- 
tally and computationally) for various flows. 

Calculations performed by this author 461 47 and other 
works 43,48 have at times shown an extreme sensitiv- 
ity to values assumed for these parameters. An exam- 
ple is taken from Ref. 47 involving calculations per- 
formed for a direct connect scramjet combustor (see Fig. ) 
tested at Wright-Patterson Air Force Base (AFRL/PRA). 
Figures and show mass-flux weighted flow properties 
through the combustor at flight conditions that corre- 
spond to Mach 4.0 and Mach 6.5 operation. Results are 
shown for various Pr, values with Sc, fixed at 0.5, and 
for several Sc, values with Pr, fixed at 0.89. The range 
of values considered is within the range of values given 
in Table 1 . As one would expect, reducing the turbulent 
Schmidt number consistently intensified combustion due 
to enhanced species diffusion processes. At the Mach 4.0 
condition, the reduction of Sc, from 0.5 to 0.25 enhanced 
turbulent mass transfer (and subsequent heat release) to 
levels that the isolator was not able to withstand, resulting 
in a potentially catastrophic un-start condition. A mod- 
est increase of Sc, from 0.5 to 0.75 reduced the turbu- 
lent mass transfer to levels that were not able to sustain 
combustion. Hence, a variation in Sc, from 0.25 to 0.75 
yielded results that covered the entire spectrum of oper- 
ability for the engine at the Mach 4.0 flight condition. 
A reduction of the turbulent Prandtl number enhanced 
the combustion process only at the higher Mach number 
state. At the Mach 4.0 condition, the heightened thermal 
diffusion processes allowed heat to be transferred away 
from the flameholding (recirculation) zones at a rate that 
was not sustainable, causing flame blow-out. These re- 
sults clearly suggest that extreme care should be taken 
when attempting to characterize these high-speed propul- 
sion devices (that contain a variety of different mix- 
ing mechanisms) with constant turbulent transport coeffi- 
cients. 


Table 1 : Turbulent Prandtl & Schmidt Number Values 


Row Field 

Pr, 

Sc, 

Planar Jets 35-38 

0.2 - 3.0 

0.1 -2.2 

Round Jets 39-41 

0.7 - 2.0 

0.1 -2.0 

Backward Facing Step 42 

0.7 - 3.0 

NA 

Jet into Cross Row 43, 44 

NA 

0.1 -0.5 

Injection Behind a Bluff Body 45 

NA 

0.2 - 0.7 


The physical mechanisms that directly influence the 
Reynolds heat and mass flux vectors can be ascertained 
by examining the transport equations for these quanti- 
ties. The transport equations that govern the Reynolds 
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heat and mass flux vectors can be written as: 
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The above expressions show that the evolution of each 
Reynolds flux vector is governed by: (I) advection, 
(II) turbulent convection, (III) source/sink due to the 
mean velocity gradient, (IV) source/sink due to the mean 
scalar gradient, (V) source/sink due to the pressure gra- 
dient, (VI) dissipation due to viscosity, and (VII) dissi- 
pation due to scalar diffusivity. Note that the Reynolds 
mass flux vector is also directly affected by the chem- 
istry (VIII). Clearly, any attempt to collapse all of these 
physical phenomena into a single gradient diffusion ef- 
fect is questionable. In fact, the literature is filled with ev- 
idence of counter-gradient diffusion effects 49-51 (i.e. tur- 
bulent diffusion of a scalar against its mean gradient) on 
the Reynolds flux vectors, particularly in pre-mixed ap- 
plications. Counter-gradient diffusion has been attributed 
to the mean pressure gradient portion of term (V) in 
Eqs. 31 and 32. 

The number of scalar transport expressions that result 
from Eqs. 3 1 and 32 is 3 x ns. Additional supporting tur- 
bulent transport equations for variances/covariances and 
their dissipation rates are also typically required to model 
the unclosed terms on the right-hand-side of these equa- 
tions. Hence, even if suitable models were developed to 
close each of the scalar flux vector equations, the number 
of additional equations introduced would greatly exceed 
the equation count given by the first order moments of 


Eqs. 5a - 5b. This fact suggests that it would be impracti- 
cal to include a full second order closure model in any 
simulation of engineering interest. To circumvent this 
difficulty, some limiting studies have invoked equilibrium 
assumptions to reduce the differential equations to alge- 
braic relationships. Other studies have coupled the gra- 
dient diffusion hypothesis with models that allow the tur- 
bulent Prandtl and/or Schmidt number to vary spatially. 

The work of Adumitroaie 52 involved the development 
and application of a complete algebraic closure for the 
Reynolds stress tensor and scalar flux vectors. The ex- 
plicit algebraic Reynolds stress model used in this effort 
was based on the closure of Taulbee 27 and included com- 
pressibility effects. The algebraic Reynolds scalar flux 
models for temperature and species composition were de- 
rived in Ref. 52 based on similar principles. The model 
neglected scalar correlations higher than second order 
and cross-correlations between temperature and compo- 
sition were neglected (as were temperature fluctuation ef- 
fects on the reaction rates). Additional transport equa- 
tions (beyond those given by Eqs. 5a - 5b) required by the 
model include the turbulent kinetic energy and its dissi- 
pation rate, the variance of temperature and its dissipation 
rate, and the variances and covariances of the ns— 1 com- 
position variables and their dissipation rates. The end re- 
sult is that nsx(ns— 1)+4 additional transport equations 
are introduced; a value that exceeds the equation count 
for the first order moments for ns> 2. Nevertheless, en- 
couraging results were obtained for a compressible mix- 
ing layer (cold flow) and planar jet when compared with 
results obtained from a fully second order moment trans- 
port model. The author noted that high shear regions were 
problematic with the model, suggesting that the highly 
nonlinear nature of the algebraic closures could pose dif- 
ficulties in complex flows. Further numerical difficulties 
associated with the use of the algebraic Reynolds mass 
flux expressions were noted by the author when chemical 
reactions were considered. 

The development of models that allow for variable tur- 
bulent Prandtl/Schmidt numbers within the context of the 
gradient diffusion hypothesis has been pursued by several 
authors. 40,41,44 The variable turbulent Prandtl number 
models tend to involve additional transport equations for 
the temperature or enthalpy/energy variance (g"g") and 
its dissipation rate (e g ). This allows an additional (in- 
dependent) turbulent time scale to be introduced into the 
definition of the thermal eddy diffusivity: 

otf : - — = C T k Xj (33) 

where C T is a model coefficient (possibly with a near- 
wall damping function) and x T is some measure of the 
turbulent thermal time scale. The thermal time scale can 
be based purely on the scalar transport variables, or a 
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mixed time scale can be defined by introducing the tur- 
bulent time scale based on the velocity field (t), i.e. 


i 



These particular choices for the thermal time scale yield 
the following expressions for the turbulent Prandtl num- 
ber: 



or 


Pr,= 


Pi 

C T kyfx 



(35) 


Variable turbulent Schmidt number models are arrived at 
in a similar fashion by integrating evolution equations for 
some measure of the composition variance (e.g. the mix- 
ture fraction variance or the sum of all species mass frac- 
tion variances) and its dissipation rate. If these quantities 
are also denoted by g"g" and e g , then the expressions re- 
lated to the eddy diffusivity of mass are obtained from 
Eqs. 33 , 34 , and 35 by replacing Pr,, C T , and x T with 
Sc t , Cfofi and x^. 


Chemical Production Rate 

The most common species production rate closures 
used for high-speed reacting flows are based on laminar- 
chemistry, eddy break-up / dissipation, 6,7 or proba- 
bility density function (PDF) 8-10 formulations. Ap- 
proaches based on a laminar-chemistry assumption sim- 
ply ignore turbulence-chemistry interactions by evaluat- 
ing the chemical source terms based on mean flow prop- 
erties. Eddy break-up models are mixing limited models 
where the chemical time scale is assumed to be limited 
by the dissipation rate of turbulent eddies. Formulations 
based on ideas borrowed from probability theory repre- 
sent perhaps the most elegant class of models for averag- 
ing the chemical source terms. However, these formula- 
tions can be considerably more expensive to invoke. 

Let a general kinetic step be denoted as follows: 


rtj+i 

X V ml^m 


m= 1 


X v'LlC m 

m= 1 


/ = l, nr 


(36) 


— W m 




ns + 1 

•v n 

n= 1 



"1 



( 37 ) 


where k t and k, are the forward and backward reac- 

Ji °i 

tion rate coefficients of reaction (typically exponen- 
tial functions of temperature), and p„ is the density of 
species “n”. The molar concentration of the third body 
constituent in Eq. 37 is defined by the following expres- 
sion: 


P ns- : 

w. 


ns+1 


X tbe ml 


m=l 


Pm 

w m 


(38) 


where tbe ml is the third body efficiency of species “m” in 
reaction “1” provided with the kinetic model. The chem- 
ical source term based on one-way global steps (some- 
times referred to as arbitrary reaction order steps) can be 
written in the following manner: 


nr . . ns / ~ \ a n \ 

£ <»> 

i=l V ' n=l V""/ 


Here, the coefficient a nl , in general, is not equal to the 
stoichiometric coefficient of species “n” in reaction “1” as 
is the case with the law of mass action. This coefficient is 
instead determined empirically using data generated from 
measurements or from a detailed kinetic mechanism. 

The species production rates are point functions (i.e. 
functions that are defined by variables at a single spatial 
and temporal location), thus they are ideally suited for 
single point PDF closures. The source terms, given by 
Eqs. 37 and 39, are a function of temperature and com- 
position only. As a result, these terms can be averaged 
by integrating the product of the species production rates 
with the joint PDF (t?) of temperature and composition 
at each spatial location, i.e. 

= J W m (t ,P] , • • • ,pns)&(t , Pj , • ■ ■ iPns) 
dfdp v ..dpns (40) 


where v ml and v" ml are the reactant and product stoichio- 
metric coefficients for species “m” in reaction C m is 
the symbol for constituent “m” (the ns + 1 constituent rep- 
resenting the third body species), and nr is the number of 
chemical reactions considered. The expressions used for 
the chemical source terms are then generally given by the 
law of mass action or empirically derived global reaction 
rate expressions. The law of mass action applies to re- 
action models that are based on elementary kinetic steps, 
and can be written as follows: 


The integration in the above expression is taken over all 
realizable values for temperature and composition, and 
the independent variables of the PDF (t and p m ) repre- 
sent the sample space of the random variables T and p m . 
The form for the joint PDF can be assumed a priori 53-57 
or by integrating the evolution equation governing the 
PDF. 8-10 Note that assumed PDF formulations based on 
mixture fraction, 53 which are popular in low-speed ap- 
plications, are seldom used in high-speed flows. These 
approaches tend to treat the reacting system as either a 
mixed-is-bumed flame sheet or assume the mixture is in 
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chemical equilibrium. Neither infinitely fast chemistry tion statistics of low-speed and high-speed mixing layers, 

assumption is appropriate in supersonic flows due to lim- This effort showed that the mixture fraction variance ex- 
ited flow residence times. traded from the high-speed compressible mixing layer 

The computational overhead associated with invoking was significantly lower than the values extracted from 
a PDF approach varies greatly depending on the partic- the low-speed case. Moreover, a measure of the mixed 

ular formulation invoked. Assumed PDF methods typi- fluid probability and the peak mean temperature were 

cally add 10-30% overhead over a laminar-chemistry cal- both higher for the high-speed mixing layer. These ob- 

culation, provided that the integrations in Eq. 40 can be servations suggest that the concept of unmixedness may 

performed analytically or through an efficient table look- be play a smaller role in high-speed reacting flows, 

up procedure. Approaches that involve the integration of Models based on the eddy dissipation concept address 
an evolution equation governing the joint PDF are con- the turbulence closure problem by assuming that the ki- 

siderably more expensive, possibly by as much as a factor netic rate is limited by the rate of mixing (on a molecu- 

of 10 or more over their laminar-chemistry counter-parts, lar scale) between fuel and oxygen carrying eddies rather 

Due to the large dimensionality of the joint PDF, solving than on the chemical time scale. In regions of high tur- 

the evolution equation with a finite difference scheme is bulence levels, the eddy lifetime is short leading to large 

not practical. 9 10 Instead, the equation is typically simu- eddy dissipation rates and more rapid molecular mixing 

lated using a Monte Carlo scheme. The additional com- than regions of lower turbulence levels. This model is 

putational cost is then dependent on the number of rep- applicable to irreversible reactions only, and is usually 

resentative sample space ensembles used for the Monte applied to a single reaction step such as: 

Carlo simulation. 

Calculations of high-speed reacting flows that have ac- v a ~ > v p Products (41) 

counted for turbulence-chemistry interactions through the , . , . 

friT -.r; f . , . r .. n f „ where the stoichiometric coefficients (v.,v„) are related 

use of PDF formulations can be found in Refs. 46 , 58 - , ... _ , . 

, . „ , .. c , , to the stoichiometric Air to Fuel mass ratio (A/F), i.e. 

61. One observation found from each of these sources v ' ' 

is that the effect of turbulence-chemistry interactions is 
relatively minor except in the vicinity of flame ignition. v A = 

Figure compares results obtained for a supersonic ax- 
isymmetric H 2 /Air burner using laminar-chemistry, as- The chemical source term based on the eddy dissipation 
sumed PDF, and evolution PDF closure approximations, concept proposed by Magnussen and Hjertager 7 is given 
The assumed PDF model invoked a Gaussian distribution jj, e f 0 u 0W i n g relationship- 
for temperature fluctuations and a multi- variate ft distri- 
bution 63 for composition fluctuations. Reaction cross- 
correlations (RCC) 61 between temperature and composi- 
tion were neglected in the model. The evolved PDF re- 
sults were obtained by integrating an equation governing where A and B are empirical constants originally set to 
the scalar probability density function for enthalpy and 4.0 and 0.5, respectively in Ref. 7. This model is popu- 

composition via a Eulerian Monte Carlo procedure. 10 lar due to its simplicity and its dependence only on first 

The similarity observed in the results extracted from order correlations (i.e. no additional transport equations 

each turbulence-chemistry closure is an outcome that is are required). Many implementations of this model also 

contrary to what is typically expected in low-speed appli- permit the chemical time scale to be considered as a lim- 

cations. Large scale mixing within turbulent eddies tends iting rate using either Eq. 37 or 39. In this scenario, the 

to “stir” the fuel and air streams rather than mix them expression that yields the smallest magnitude for the re- 

at a scale small enough for chemical reactions to take action rate is the expression used to compute the source 

place. Hence, the average mass fractions within a control term. This additional limit discourages chemical reac- 

volume larger than the eddy would suggest that the two tions in cold regions of the flowfield. One clear advan- 

streams are well mixed, but in reality the two stream may tage of the pure model given by Eq. 43 is that reaction 

still be segregated within the eddy. This phenomenon, re- rate coefficients, which are often not available for com- 

ferred to in the literature as unmixedness, leads to large plex fuels, are not required. This model also requires a 

scalar covariance levels, and tends to substantially reduce minimal number of species transport equations, making 

the magnitude of the species production rate as compared it computationally efficient. Moreover, the use of a single 

with results based on a laminar-chemistry treatment. One time scale for reaction (t) alleviates much of the stiff- 

explanation as to why this phenomenon is not as preva- ness involved with more complex chemical systems. The 

lent in high-speed flows was described in Ref. 38. In major drawback of this model is that the details of the 

this work, LES was used to examine the scalar fluctua- chemical processes are neglected. Consequently, models 


— yy ( y" _ v ^ —MIN — B P p 

~ m V m m ) x ’ v.W. ” VpWp 
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of this type should never be used to predict lean blow-out 
limits or combustor ignition characteristics. 

CONCLUSIONS 

The use of steady-state RANS models has been and 
will continue to be the tool of choice for modeling com- 
pressible reacting flows for high-speed commercial and 
military applications. Even with expected increases in 
computer speed, the role of LES will likely be limited to 
idealized component analysis, or to scenarios where flow 
unsteadiness is of special concern. Hence, improvements 
to modeling approaches for compressible reacting flows 
within a RANS framework offers the greatest potential 
advancement to CFD practitioners. Of the issues raised 
in this document, improvements to the turbulent scalar 
transport models are likely to reap the most benefits. The 
simple gradient diffusion models with constant transport 
coefficients have proven to be particularly troublesome. 
When one considers the limited residence times associ- 
ated with scramjet engine flows (typically on the order of 
one millisecond), it is not surprising how even a small dis- 
crepancy in mixing rate prediction can lead to large devi- 
ations in combustor performance. Considerably more at- 
tention has historically been given to higher order models 
for closure of the Reynolds stress tensor (at least in low- 
speed applications). The use of models from this class 
is envisioned to improve upon predictions of hypersonic 
inlet and isolator flows which are dominated by shock- 
induced separation and Reynolds stress anisotropies. Lin- 
ear Reynolds stress models are certainly not capable of 
predicting the latter of these flow scenarios. Calculations 
to date have yet to show a first order need for the ad- 
vancement of turbulence-chemistry interaction models in 
high-speed applications; although this issue has not yet 
been thoroughly addressed. This observation is in stark 
contrast to low-speed reacting flows where these models 
are required to avoid a significant over-prediction of the 
mean temperature field. At this time, imprecise results 
associated with the modeling of turbulence-chemistry in- 
teractions tend to be overshadowed by inaccuracies in tur- 
bulent scalar transport predictions. 
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Figure 1: Measured (left) and computed ( k-a> with Boussinesq) pitot pressure profiles at x/h=40 



Figure 2: Measured (left) and computed (explicit algebraic Reynolds stress) pitot pressure profiles at x/h=40 



Figure 3: Measured and computed wall pressure distributions along the center of the duct 
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Figure 7: Mean H 2 0 mole fraction comparisons with measurements (Cheng et al .) obtained from laminar chemistry, 
assumed PDF, and evolved PDF models 
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